The present disclosure relates generally to imaging systems and particularly to systems and methods of reconstructing an image using iterative techniques.
Traditionally, images have been reconstructed from Computed Tomography (CT) data using so-called direct reconstruction algorithms such as filtered back projection (FBP) or convolution back projection (CBP). Recently, iterative reconstruction (IR) algorithms have received renewed attention for the reconstruction of CT images. The major advantage of IR algorithms is that they can more accurately model the measurements obtained from real CT systems. This is particularly true for helical CT systems with multi-slice detectors because these systems produce projection measurements that pass obliquely through the 2-D reconstructed image planes. By more accurately modeling these projections, IR algorithms can produce reconstructions with higher quality, lower noise, and fewer artifacts.
For example, consider a helical scan CT system. The 3-D volume to be reconstructed can be represented by an array of N discrete voxels xi where i is the index of the voxel 3-D position. The value xi may specify the unknown density of the voxel. Furthermore, let x=[x1, x2, . . . xN] be a vector containing the unknown density of each voxel in the reconstruction. In this case, x represents the full 3-D reconstruction volume. During the CT scanning process, projection data are measured for M different projection lines through the object. The different projection lines are typically measured for a wide variety of positions and angles through the object. The value of the jth projection through the object is denoted by yj and the vector of all measurements is denoted byy=[y1, y2, . . . , yM].
The objective of IR algorithms is to determine the unknown value of x by searching for the value of the vector x that best matches the measured data. Typically, this is done by minimizing a cost function of the form:
                              x          ^                =                  arg          ⁢                                          ⁢                                    min              x                        ⁢                          {                                                ∑                                      m                    =                    0                                    M                                ⁢                                                      w                    m                                    ⁢                                                                                                                                    y                          m                                                -                                                                              F                            m                                                    ⁡                                                      (                            x                            )                                                                                                                                      2                                                              }                                                          equation        ⁢                  -                ⁢        1            where {circumflex over (x)} is the value of the variable x which achieves the minimum of the function and wm is a weighting factor, which can be based upon statistical properties of the algorithm. This cost function can be minimized in a variety of manners using optimization methods such as iterative coordinate descent, expectation maximization, conjugate gradient, a gradient-based method, or any number of alternative techniques.
In practice, the solution to equation-1 often includes greater levels of noise than desired. This noisiness may result when there are too few measurements, when the quality of the measurements is poor, or When the available projection angles and locations do not give sufficient information about x to properly reconstruct it. This problem can be addressed by adding an additional stabilizing function S(x) to the cost function being minimized. This results in the regularized inverse
                              x          ^                =                  arg          ⁢                                    min              x                        ⁢                          {                                                                    ∑                                          m                      =                      0                                        M                                    ⁢                                                            w                      m                                        ⁢                                                                                                                                                y                            m                                                    -                                                                                    F                              m                                                        ⁡                                                          (                              x                              )                                                                                                                                                  2                                                                      +                                  S                  ⁡                                      (                    x                    )                                                              }                                                          equation        ⁢                  -                ⁢        2            
Typically, the function S(x) is chosen to be a quadratic function with the formS(x)=xtHxwhere H is a symmetric and positive definite or positive semi-definite matrix. The particular choice of the functional S(x) can have a substantial effect on the quality of reconstructions produced by IR algorithms.
A major challenge of IR is the computation time and computational resources required to complete a reconstruction. Because IR has been studied for other types of reconstruction problems, a variety of methods that have been proposed for computing the solution to equation (1). Some of these methods include ordered subset expectation maximization (OSEM), preconditioned conjugate gradient (PCG), and iterative coordinate descent (ICD). All of these methods perform the minimization required in equation (1) by iteratively computing a forward projection Fm(x), determining the error between the forward projections and the measured projections ym, adjusting the values in the image x, and then repeating this process until the result is sufficiently close to the desired minimum.
Of the available IR reconstruction methods, ICD has been found to provide convergence rapidly, or with few required iterations. The ICD algorithm differs from the others in that it changes or updates a single voxel at a time, and sequentially visits all voxels in the reconstruction volume. An N-dimensional optimization problem is reduced to a sequence of one-dimensional greedy optimizations, where at each step the largest amount of change in a voxel which minimizes the 1D cost function is selected as the step to descend the global objective function. In the conventional ICD algorithm, every voxel is updated exactly once for each iteration of the algorithm. The order of voxel selection for 1D updates is important to reduce the correlation between successive updates and take as large a step to descend the global cost function as possible with each voxel update. For instance, it is generally known in the art that randomized voxel selection for ICD updates results in faster convergence than simple raster scan ordering.
While the convergence of ICD is rapid in terms of the number of required iterations, the computational requirements of 3D ICD for helical multi-slice CT are still very demanding. Various methods have been investigated for reducing the computational requirements to ICD make faster. One method is to focus computations upon a particular region-of-interest (ROI). Another method is to perform reconstructions at different resolutions, typically going from a coarse resolution representation of the reconstructed image to finer resolutions reconstructions. While these methods are effective in certain applications, they have disadvantages. ROI reconstruction only produces a reconstruction in a limited spatial region. If it is necessary to view the reconstruction in a larger area, then ROI reconstruction is unsatisfactory. Multi-resolution algorithms can speed reconstruction however they lead to complex algorithms, which can be difficult to implement and maintain.
Accordingly, there is a need in the art for an iterative reconstruction arrangement that overcomes these drawbacks.